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We generalize the Kuramoto model of coupled oscillators 
to allow time-delayed interactions. New phenomena include 
bistability between synchronized and incoherent states, and 
unsteady solutions with time-dependent order parameters. 
We derive exact formulas for the stability boundaries of the 
incoherent and synchronized states, as a function of the de- 
lay, in the special case where the oscillators are identical. The 
experimental implications of the model are discussed for pop- 
ulations of chirping crickets, where the finite speed of sound 
causes communication delays, and for physical systems such 
as coupled phase-locked loops or lasers. 
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The Kuramoto model of coupled oscillators is one of 
the most celebrated systems in nonlinear dynamics. It 
was originally developed as an analytically tractable ver- 
sion of Winfree's mean-field model for large populations 
of biological oscillators 0|, such as groups of chorusing 
crickets ||], flashing fireflies ||, and cardiac pacemaker 
cells [|. In a beautiful analysis, Kuramoto showed that 
the model exhibits a spontaneous transition from inco- 
herence to collective synchronization, as the coupling 
strength is increased past a certain threshold ||. The 
model has since been analyzed more deeply and extended 
in various ways |6|-|To| . It has also been linked to several 
physical problems, including Landau damping in plas- 
mas l|, the dynamics of Josephson junction arrays [0 , 
bubbly fluids [Q , and coupled Brownian ratchets |13|]. 

Here we explore the effects of time delay on the dynam- 
ics of the Kuramoto model. In the past, delay has often 
been neglected in models of coupled oscillators. In many 
cases this approximation is physically justified, and in 
all cases it simplifies the mathematics. But recently sev- 
eral authors have begun to investigate oscillator systems 
where delays are not negligible fl^ , |l5| , motivated by neu- 
ral networks where synaptic, dendritic, and propagation 
delays can be significant. Other authors have considered 
delays in systems of limit-cycle oscillators |lq] , with ap- 
plications to arrays of lasers and microwave oscillators. 

Intuitively, the problem is similar to that faced by the 
fans sitting in an enormous football stadium, all of whom 
(we suppose) are trying to clap in unison. Even if every- 
one were successfully clapping in perfect synchrony, it 
would not sound that way to the fans themselves, as the 
applause coming from far across the field would be sig- 
nificantly delayed, because of the finite speed of sound. 



Nevertheless, we show that perfect synchrony is possi- 
ble in the Kuramoto model with time delay, if all oscilla- 
tors are identical. In fact, there can be several different 
synchronized states, and they can co-exist with a stable 
incoherent state where the oscillators are completely dis- 
organized. These multistabilities are qualitatively new: 
they do not occur in the original Kuramoto model. 

We consider a system of phase oscillators with noisy, 
randomly distributed intrinsic frequencies, and with de- 
layed mean-field coupling: 
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for i = 1, . . . ,N. Here 9i(t) is the phase of the ith oscilla- 
tor at time t, and u)i is its intrinsic frequency, randomly 
drawn from a probability density g(w) with mean loq- The 
white noise £i(f) represents frequency fluctuations at an 
effective temperature D > 0, and is defined by the ensem- 
ble averages < >= 0, < &(»)£/'(*) >= 2D5 lj 5(s-t). 
In the global coupling term, K > is the coupling 
strength, r > is the delay, and a is a phase frustra- 
tion parameter. This model reduces to the Kuramoto 
model H if r = 0, a = 0, and D = 0, and to the mean- 
field XY model if r = 0, a = 0, and the oscillators are 
identical, i.e., g(w) = 5(lu — ujq). For r = 0, the separate 
effects of frustration a and noise D have been studied by 
Sakaguchi and Kuramoto ||. 

As the one-parameter family of rotating-frame trans- 
formations 9i(t) — > 6i(t) — fit, u>i — ► u>i — Q, a — > a + 
leave Eq. (|l|) invariant for any fl, we may assume a = 
without loss of generality — except if r = 0, which we for- 
bid. (This restriction is merely for convenience. All our 
results are well-behaved as r — > and converge to those 
obtained by setting r = from the start.) Moreover, 
since Eq. ([!]) is invariant under the reflection u>i — > — a;,-, 
8i — > — 6i, a — > — a. it suffices to consider ujq > 0. 

It is often helpful to describe the macroscopic state 
of the system in terms of the complex order parameter 
fl(f) e #W = j. J2f =1 e^ (t ) introduced by Kuramoto [§. 
Here R(t) measures the system's phase coherence. In par- 
ticular, R — 1 if all the oscillators are in phase, whereas 
R = if the oscillators are scattered around the unit 
circle with their centroid at the origin. 

Our first analytical result concerns the stability of the 
incoherent state for the infinite- N limit of Eq. ([!]). We 
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rewrite the model as a Fokker-Planck equation for the 
density p(0,u),t) of oscillators currently at phase 0, with 
intrinsic frequency u>. Because the method is standard 
|j],D, we omit the details. The only new feature here 
is that the drift velocity field inherits the time delay 
in Eq. (Q). When the Fokker-Planck equation is lin- 
earized about the incoherent state (the stationary density 
p = 1/2tt), we find Jl7| that its continuous spectrum is 
{— D — iu) | uj G supp(g)}. Hence for D > 0, the continu- 
ous spectrum corresponds to damped modes and there- 
fore the stability of the incoherent state is determined 
solely by the discrete eigenvalues. But when D = 0, the 
continuous spectrum is pure imaginary and corresponds 
to neutrally stable rotating waves in the full system. In 
this case, the incoherent state can never be linearly sta- 
ble: it is either unstable or neutral, depending on the 
discrete eigenvalues. These eigenvalues A satisfy |17[] 
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This implicit formula for A is exact but difficult to 
analyze for arbitrary g(uj), so we consider the case of 
identical oscillators to gain some insight. Even this case 
turns out to be far from trivial. If g(u) = S(lu — ujq), Eq. 
(0) can be simplified to the transcendental equation 



(p + ir — z)e z + q = 0, 



(3) 



where p = —Dt < 0,r = — ujqt, q = Kt/2, and z = Ar. 
Then the stability of the incoherent state depends on 
whether all roots of Eq. (||) satisfy Re(z) < (in which 
case we will say "all eigenvalues are stable" , for brevity) . 

By the transformations z — > z + imr, q — > (— l) n q, r — > 
r + nit, and z — > z*, r — > — r, we may assume r G [0, 7r/2] 
in Eq. (|). For r = 0, Hayes proved @|| that all 
eigenvalues are stable if and only if p < 1 and p < —q < 
\fp 2 + yi 2 , where y\ is the unique zero of p sin y — y cos y 
in (0,7r). Using results of Pontryagin as in Ref. Q, we 
can show R|o| that for r G (0,7r/2], 



• if p = 0, then all eigenvalues are stable if and only 
if r - n/2 < q < 0. 

• if p < 0, then all eigenvalues are stable if and only 
if - y/p 2 + (y 2 — r) 2 < q < \/ p 2 + (yi - r) 2 , where 
?/i and y 2 are the unique zeroes of p sin y + (r — 
y)cosy in (0,r) and (tt/2, it) respectively. 

These conditions are exact but still opaque, so we sim- 
plify the model further for illustration. Suppose there is 
no noise (D = 0). Then we find JItJ that the incoherent 
state is neutrally stable precisely when 
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with m being an arbitrary positive integer. 



Figure [j] shows that this analytical result agrees with 
numerical simulations, even for as few as N — 12 oscil- 
lators. Although the incoherent state is neutrally stable 
in the grey region, we observe numerically that R.(t) — » 
exponentially fast, as in Landau damping ||. In this 
sense, incoherence is stable in the grey region. 
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FIG. 1. Stability region for the incoherent state, with 
g(uj) = S(u - tJo), luo = k/2 D = 0,N = 12. Black 
curves: theoretical boundaries (0) for the infinite-iV limit; 
Grey area: results from numerical integration using a sixth 
order Adams-Bashforth-Moulton scheme, with fixed stepsize 
dt = t/20, and with the corrector formula iterated for con- 
vergence and stability. Initially, all the phases were evenly 
spaced, and the symmetry was broken by adding O(10~ 10 ) 
random perturbations. The incoherent state was judged as 
unstable if R(t) > 10 -7 at a final time of t = 800t. 

Continuing with the instructive case of noiseless, iden- 
tical oscillators, we now consider the possibility of perfect 
synchrony: 0j(t) = 0(t) for all i. We restrict our attention 
to a particular class of such solutions, namely uniform ro- 
tations: 8(t) = fit + j3. Self-consistency then requires 

fl = uj - K sin(fir) (5) 

for such solutions to exist, and linearization fl7|| imposes 

cos(r2r) > (6) 

as the condition for their orbital stability. 

If we graph both sides of Eq. (JsJ) as functions of f2, we 
see that for all sufficiently large K, there exist multiple 
stable synchronized states |2lJ , as Eq. (|^) has non- unique 
solutions satisfying (||). We can also see that stable syn- 
chrony is impossible for certain combinations of r and K. 
The problem reduces to characterizing the two-parameter 
family of lines of negative slope that intersect the sine 
function on its descending limbs. We find that stable 
synchronized states do not exist if and only if 
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with m being an arbitrary positive integer. 
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These zones of forbidden synchrony are shown in black 
in Fig. ||. For comparison, they are overlaid on top of the 
earlier grey regions (Fig. ||) where incoherence is stable. 
The black regions fit neatly inside the grey; they have the 
same base and half the height. The exposed parts of the 
grey regions correspond to bistability: stable incoherence 
coexists with stable synchrony, and hysteresis can occur. 
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FIG. 2. Stability regions of the incoherent state ((^])) and 
the synchronized states ((Q)) for g(co) = 8(u — c^o), D = 0. 
White region: one or more stable synchronized states exist 
but the incoherent state is unstable; Black region: incoherence 
is stable but synchrony is not; Grey region: one or more stable 
synchronized states coexist with stable incoherence. 

Numerical simulations reveal windows in the bistable 
regions where R can be time-periodic (Fig. |^) . In the ex- 
ample shown, period doubling occurs as K increases, but 
seems to be truncated beyond period 16 (not shown); 
after that, a synchronized state apparently takes over, 
suppressing further period doubling (Fig. ||(d)). Such 
unsteady behavior is a consequence of the delay; in the 
standard Kuramoto model, numerical experiments show 
that R(t) always approaches a constant value if g(u>) is 
unimodal and symmetric (although this has never been 
proven) . Oscillator configurations with two or more clus- 
ters |22| cause the unsteady behavior seen here. All the 
clusters move with the same average velocity, but their 
separation is periodically modulated. 

So far we have concentrated on identical oscillators. 
To check how the results would be modified for other fre- 
quency distributions, we have considered the Lorcntzian 
distribution g(iv) = (7/71) (j 2 + (u — cj ) 2 ) -1 - Then by a 
remarkable coincidence [jlTJ , Eq. (||) can again be reduced 
to Eq. (||), but now with p = — (7 + D)t, r = —ujqt, q = 
Kt/2, and z — At. The critical coupling is given by 

K c = 2{j + D)sec(tt c T), (8) 

where fl c = luq — (7 + D) tan(Sl c r). Figure [| plots the 
corresponding region where incoherence is stable. It re- 
sembles Fig. 0, with a series of evenly spaced peaks (at 
ujqt = (2n + l)7r) that decrease in height. The main dif- 
ference is that the distributed frequencies produce some 
rounding of the boundary, and lift it off the r-axis so that 
it now has minima 2(7 + D) at luqt = 2mt. 
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FIG. 3. Time series showing nonsteady order parameter 

R(t), with g(u) = S(lu - uj ), uj = tt/2, D = 0, r = 2, TV = 24. 

(a) K = 1.3 : period-2 oscillation; (b) K — 1.4 : period-4 

oscillation; (c) K = 1.44 : period-8 oscillation; (d) K = 1.475 : 

R(t) — > 1 after a periodic transient. 



The bistability found earlier also has a counterpart in 
the Lorentzian case (Fig. ||](a)). Partially locked states 
(which may not be unique) replace the earlier in-phase 
states, but otherwise the story is unchanged ^3|. Thus, 
the case of identical oscillators captures the essential fea- 
tures introduced by delay. 

Our final result concerns the bifurcation at K = K c , 
where the incoherent state becomes unstable. We have 
adapted the two-timing method of Ref. K ] to handle the 
delay-differential equations |l]) . We find [ 17j that gener- 
ically, for D > and arbitrary g(u>) 7 a Hopf bifurcation 
occurs at K c [^ij , giving rise to a partially locked state, or 
in the density description, a rotating wave with constant 
coherence R = 0(-J\K — K c \). This bifurcation may be 
subcritical (Fig. [|(a)) or supercritical (Fig. ||(b)). 

Experimental tests of the model may be possible in 
arrays of phase- locked loop s p5| , relativistic magnetrons 
p6[ , or solid-state lasers [p7[, as they are all approxi- 
mately governed by coupled Adler equations [^8| similar 
in form to Eq. (0). The delay and the coupling strength 
are both natural control parameters, and perhaps one 
could try to map out the stability boundaries, look for 
hysteresis between incoherence and synchrony, etc. Our 
model may also help to explain how crickets can synchro- 
nize their chirps , despite the time delays caused by the 
speed of sound. Crickets listen to each other's chirps and 
adjust their own timing according to a phase response 
curve Q| . The propagation delay between two crickets 3 
meters apart is about 10 msec. This is short compared 
to the chirp period (300-500 msec). Our results suggest 
that delay effects become significant only in the first peak 
in Fig. |2[ i.e., for delays near half the period of the oscil- 
lation. Thus the delays that crickets actually encounter 
in the field are probably negligible as far as synchrony is 
concerned. It would be interesting to try lab experiments 
on crickets interacting via chirp signals whose delay and 
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FIG. 4. Stability region of the incoherent state for 
Lorentzian g(u>) with ujo — 3, 7 = 0.1, D = 0, N = 3600. 
Black curves: theoretical boundary Eq. (g); Grey strips: nu- 
merical results using the same method as in Fig. oL Insets: (a) 
Subcritical Hopf bifurcation of the incoherent state at r = 1. 
Arrows indicate a hysteresis loop between stable incoherence 
and a stable partially locked state with R close to 1. (b) Su- 
percritical Hopf bifurcation of the incoherent state at r = 2. 
A stable partially locked state grows continuously from the 
incoherent state with R — 0(yK — K c ). 



amplitude can be electronically manipulated. 
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